
#define PARALLEL
// #define SSE2
#define USE_DOUBLE_PRECISION
#define NO_RANDOM_PER_MDP_SITE
#include "fermiqcd.h"
#include "fermiqcd_lanczos.h"

int main(int argc, char** argv) {
  mpi.open_wormholes(argc,argv);
  define_base_matrices();
  int i, alpha, iter, box[]={20,10,10,10};
  mdp_lattice lattice(4,box);
  site        x(lattice);
  gauge_field U(lattice,3);
  fermi_field psi(lattice,3);
  Complex ab;

  if(argc < 3) error("Usage: test_lanczos.exe [n] [gauge] [trifile]");
  
  if(ME==0) iter=(int) val(argv[1]);
  mpi.broadcast(iter,0);
  
  U.load_as_double(argv[2]);
  

  float p;
  forallsites(x) {
    for(i=0; i<4; i++) {
      p=max(hermitian(U(x,i))-inv(U(x,i)));
      if(p>1.0e-5) {
	printf("ME=%i x=(%i %i %i %i) mu=%i det=%f precision=%e\n",
	       ME, x(0), x(1), x(2), x(3), i, abs(det(U(x,i))), p);
      }
    }
  }

  if(ME==0) printf("\nComputing average plaquette\n");
  mdp_real plaquette=average_plaquette(U);
  if(ME==0) printf("done. Average_plaquette = %e\n", plaquette);

  psi=0;
  if(on_which_process(lattice, 0,0,0,0)==ME) {
    x.set(0);
    for(alpha=0; alpha<4; alpha++)
      for(i=0; i<3; i++)
	psi(x,alpha,i)=1;
  }

  coefficients coeff;
  coeff["kappa_s"]=coeff["kappa_t"]=0.1415;
  coeff["c_{sw}"]=1.57;

  if(ME==0) printf("\nComputing Chromo-Electro-Magnetic field...\n");
  compute_em_field(U);
  if(ME==0) printf("done.\n");


#if defined(SSE2)
  default_fermi_action=FermiCloverActionSSE2::mul_Q;
#else
  default_fermi_action=FermiCloverActionFast::mul_Q;
#endif
  
  double time=mpi.time();
  FILE *fp=0;
  if(ME==0 && argc==4) fp=fopen(argv[3], "w");
  for(i=0; i<iter; i++) {
    ab=Lanczos<fermi_field>::step(psi, U, coeff);
    if(ME==0 && argc==4) {
      if(i==0) printf("\nStep\tAlpha\t\tBeta\t\tTime\n");
      printf("%i\t%+20.14e\t%+20.14e\t%e\n", 
	     i, real(ab), imag(ab), mpi.time()-time);
      if(fp) fprintf(fp,"%i\t%+20.14e\t%+20.14e\t%e\n", 
		     i, real(ab), imag(ab), mpi.time()-time);
    }
  }
  if(ME==0 && fp) fclose(fp);    

  time=mpi.time()-time;

  if(ME==0)
    printf("\nTime per site per sweep = %e\n",
	   time/(i*lattice.size()));

  if(ME==0)
    printf("\nTotal time = %e\n", time);

  if(ME==0)
    printf("%f %f %f %f %f\n",
	   t1+t2+t3+t4, 
	   t1/(t1+t2+t3+t4), t2/(t1+t2+t3+t4),t3/(t1+t2+t3+t4),t4/(t1+t2+t3+t4));
  
  mpi.close_wormholes();
  return 0;
}

